{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W0D5_Statistics/W0D5_Tutorial1.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# Neuromatch Academy: Precourse Week, Day 5, Tutorial 1\n",
    "# Introduction to Probability and Statistics\n",
    "\n",
    "__Content creators:__ Ulrik Beierholm\n",
    "\n",
    "__Content reviewers:__ Ethan Cheng, Manisha Sinha\n",
    "\n",
    "Name Surname, Name Surname. This includes both reviewers and editors. Add reviewers first then editors (paper-like seniority :) ).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "#Tutorial Objectives\n",
    "\n",
    "We will here cover the basic ideas from probability and statistics, as a reminder of what you have hopefully previously been taught. These ideas will be important for almost every one of the following topics covered in the course. \n",
    "\n",
    "There are many additional topics within probability and statistics that we will not cover as they are not central to the main course. We also do not have time to get into a lot of details, but hopefully this will help you recall material you have previously encountered.\n",
    "\n",
    "\n",
    "By completing the exercises in this tutorial, you should:\n",
    "* get some intuition about how stochastic randomly generated data can be\n",
    "* understand how to model data using simple probability distributions\n",
    "* understand the difference between discrete and continuous probability distributions\n",
    "* be able to plot a Gaussian distribution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Setup\n",
    "Make sure to run this before you get started."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:25.743471Z",
     "iopub.status.busy": "2021-06-02T20:55:25.742790Z",
     "iopub.status.idle": "2021-06-02T20:55:27.314357Z",
     "shell.execute_reply": "2021-06-02T20:55:27.313138Z"
    }
   },
   "outputs": [],
   "source": [
    "# Imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "from scipy.stats import norm  # the normal probability distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:27.321648Z",
     "iopub.status.busy": "2021-06-02T20:55:27.321064Z",
     "iopub.status.idle": "2021-06-02T20:55:27.516995Z",
     "shell.execute_reply": "2021-06-02T20:55:27.516227Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Figure settings\n",
    "import ipywidgets as widgets    # interactive display\n",
    "from ipywidgets import interact, fixed, HBox, Layout, VBox, interactive, Label, interact_manual\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")\n",
    "#plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/NMA2020/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:27.532321Z",
     "iopub.status.busy": "2021-06-02T20:55:27.520245Z",
     "iopub.status.idle": "2021-06-02T20:55:27.534496Z",
     "shell.execute_reply": "2021-06-02T20:55:27.533920Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Plotting Functions\n",
    "\n",
    "def plot_random_sample(x, y, figtitle = None):\n",
    "  \"\"\" Plot the random sample between 0 and 1 for both the x and y axes.\n",
    "\n",
    "    Args:\n",
    "      x (ndarray): array of x coordinate values across the random sample\n",
    "      y (ndarray): array of y coordinate values across the random sample\n",
    "      figtitle (str): title of histogram plot (default is no title)\n",
    "\n",
    "    Returns:\n",
    "      Nothing.\n",
    "  \"\"\"\n",
    "  fig, ax = plt.subplots()\n",
    "  ax.set_xlabel('x')\n",
    "  ax.set_ylabel('y')\n",
    "  plt.xlim([-0.25, 1.25]) # set x and y axis range to be a bit less than 0 and greater than 1\n",
    "  plt.ylim([-0.25, 1.25])\n",
    "  plt.scatter(dataX, dataY)\n",
    "  if figtitle is not None:\n",
    "    fig.suptitle(figtitle, size=16)\n",
    "  plt.show()\n",
    "\n",
    "def plot_random_walk(x, y, figtitle = None):\n",
    "  \"\"\" Plots the random walk within the range 0 to 1 for both the x and y axes.\n",
    "\n",
    "    Args:\n",
    "      x (ndarray): array of steps in x direction\n",
    "      y (ndarray): array of steps in y direction\n",
    "      figtitle (str): title of histogram plot (default is no title)\n",
    "\n",
    "    Returns:\n",
    "      Nothing.\n",
    "  \"\"\"\n",
    "  fig, ax = plt.subplots()\n",
    "  plt.plot(x,y,'b-o', alpha = 0.5)\n",
    "  plt.xlim(-0.1,1.1)\n",
    "  plt.ylim(-0.1,1.1)\n",
    "  ax.set_xlabel('x location')\n",
    "  ax.set_ylabel('y location')\n",
    "  plt.plot(x[0], y[0], 'go')\n",
    "  plt.plot(x[-1], y[-1], 'ro')\n",
    "\n",
    "  if figtitle is not None:\n",
    "    fig.suptitle(figtitle, size=16)\n",
    "  plt.show()\n",
    "\n",
    "def plot_hist(data, xlabel, figtitle = None, num_bins = None):\n",
    "  \"\"\" Plot the given data as a histogram.\n",
    "\n",
    "    Args:\n",
    "      data (ndarray): array with data to plot as histogram\n",
    "      xlabel (str): label of x-axis\n",
    "      figtitle (str): title of histogram plot (default is no title)\n",
    "      num_bins (int): number of bins for histogram (default is 10)\n",
    "\n",
    "    Returns:\n",
    "      count (ndarray): number of samples in each histogram bin\n",
    "      bins (ndarray): center of each histogram bin\n",
    "  \"\"\"\n",
    "  fig, ax = plt.subplots()\n",
    "  ax.set_xlabel(xlabel)\n",
    "  ax.set_ylabel('Count')\n",
    "  if num_bins is not None:\n",
    "    count, bins, _ = plt.hist(data, bins = num_bins)\n",
    "  else:\n",
    "    count, bins, _ = plt.hist(data, bins = np.arange(np.min(data)-.5, np.max(data)+.6)) # 10 bins default\n",
    "  if figtitle is not None:\n",
    "    fig.suptitle(figtitle, size=16)\n",
    "  plt.show()\n",
    "  return count, bins\n",
    "\n",
    "def my_plot_single(x, px):\n",
    "  \"\"\"\n",
    "  Plots normalized Gaussian distribution\n",
    "\n",
    "    Args:\n",
    "        x (numpy array of floats):     points at which the likelihood has been evaluated\n",
    "        px (numpy array of floats):    normalized probabilities for prior evaluated at each `x`\n",
    "\n",
    "    Returns:\n",
    "        Nothing.\n",
    "  \"\"\"\n",
    "  if px is None:\n",
    "      px = np.zeros_like(x)\n",
    "\n",
    "  fig, ax = plt.subplots()\n",
    "  ax.plot(x, px, '-', color='C2', LineWidth=2, label='Prior')\n",
    "  ax.legend()\n",
    "  ax.set_ylabel('Probability')\n",
    "  ax.set_xlabel('Orientation (Degrees)')\n",
    "\n",
    "def plot_gaussian_samples_true(samples, xspace, mu, sigma, xlabel, ylabel):\n",
    "  \"\"\" Plot a histogram of the data samples on the same plot as the gaussian\n",
    "  distribution specified by the give mu and sigma values.\n",
    "\n",
    "    Args:\n",
    "      samples (ndarray): data samples for gaussian distribution\n",
    "      xspace (ndarray): x values to sample from normal distribution\n",
    "      mu (scalar): mean parameter of normal distribution\n",
    "      sigma (scalar): variance parameter of normal distribution\n",
    "      xlabel (str): the label of the x-axis of the histogram\n",
    "      ylabel (str): the label of the y-axis of the histogram\n",
    "\n",
    "    Returns:\n",
    "      Nothing.\n",
    "  \"\"\"\n",
    "  fig, ax = plt.subplots()\n",
    "  ax.set_xlabel(xlabel)\n",
    "  ax.set_ylabel(ylabel)\n",
    "  # num_samples = samples.shape[0]\n",
    "\n",
    "  count, bins, _ = plt.hist(samples, density=True)\n",
    "\n",
    "  plt.plot(xspace, norm.pdf(xspace, mu, sigma),'r-')\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "# Section 1: Stochasticity and randomness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:27.542550Z",
     "iopub.status.busy": "2021-06-02T20:55:27.541965Z",
     "iopub.status.idle": "2021-06-02T20:55:27.645913Z",
     "shell.execute_reply": "2021-06-02T20:55:27.645410Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 1: Stochastic world\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"UCzDNJvChv8\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 1.1: Intro to Randomness\n",
    "\n",
    "Before trying out different probability distributions, let's start with the simple uniform distribution, U(a,b), which assigns equal probability to any value between a and b.\n",
    "\n",
    "To show that we are drawing a random number $x$ from a uniform distribution with lower and upper bounds $a$ and $b$ we will use this notation:\n",
    "$x \\sim U(a,b)$. Alternatively, we can say that all the potential values of $x$ are distributed as a uniform distribution between $a$ and $b$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coding Exercise 1.1: Create randomness\n",
    "\n",
    "Numpy has many functions and capabilities related to randomness.  We can draw random numbers from various probability distributions. For example, to draw 5 uniform numbers between 0 and 100, you would use `np.random.uniform(0, 100, size = (5,))`. \n",
    "\n",
    " We will use `np.random.seed` to set a specific seed for the random number generator. For example, `np.random.seed(0)` sets the seed as 0. By including this, we are actually making the random numbers reproducible, which may seem odd at first. Basically if we do the below code without that 0, we would get different random numbers every time we run it. By setting the seed to 0, we ensure we will get the same random numbers. There are lots of reasons we may want randomness to be reproducible. In NMA-world, it's so your plots will match the solution plots exactly! \n",
    "\n",
    "```\n",
    "np.random.seed(0)\n",
    "random_nums = np.random.uniform(0, 100, size = (5,))\n",
    "```\n",
    "\n",
    "Below, you will complete a function `generate_random_sample` that randomly generates `num_points` $x$ and $y$ coordinate values, all within the range 0 to 1. You will then generate 10 points and visualize.\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:27.652758Z",
     "iopub.status.busy": "2021-06-02T20:55:27.651820Z",
     "iopub.status.idle": "2021-06-02T20:55:27.754798Z",
     "shell.execute_reply": "2021-06-02T20:55:27.754181Z"
    }
   },
   "outputs": [],
   "source": [
    "def generate_random_sample(num_points):\n",
    "  \"\"\" Generate a random sample containing a desired number of points (num_points)\n",
    "  in the range [0, 1] using a random number generator object.\n",
    "\n",
    "  Args:\n",
    "    num_points (int): number of points desired in random sample\n",
    "\n",
    "  Returns:\n",
    "    dataX, dataY (ndarray, ndarray): arrays of size (num_points,) containing x\n",
    "    and y coordinates of sampled points\n",
    "\n",
    "  \"\"\"\n",
    "\n",
    "  ###################################################################\n",
    "  ## TODO for students: Draw the uniform numbers\n",
    "  ## Fill out the following then remove\n",
    "  raise NotImplementedError(\"Student exercise: need to complete generate_random_sample\")\n",
    "  ###################################################################\n",
    "\n",
    "  # Generate desired number of points uniformly between 0 and 1 (using uniform) for\n",
    "  #     both x and y\n",
    "  dataX = ...\n",
    "  dataY = ...\n",
    "\n",
    "  return dataX, dataY\n",
    "\n",
    "# Set a seed\n",
    "np.random.seed(0)\n",
    "\n",
    "# Set number of points to draw\n",
    "num_points = 10\n",
    "\n",
    "# Draw random points\n",
    "dataX, dataY = generate_random_sample(num_points)\n",
    "\n",
    "# Visualize\n",
    "plot_random_sample(dataX, dataY, \"Random sample of 10 points\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:27.761971Z",
     "iopub.status.busy": "2021-06-02T20:55:27.760047Z",
     "iopub.status.idle": "2021-06-02T20:55:28.041283Z",
     "shell.execute_reply": "2021-06-02T20:55:28.041765Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def generate_random_sample(num_points):\n",
    "  \"\"\" Generate a random sample containing a desired number of points (num_points)\n",
    "  in the range [0, 1] using a random number generator object.\n",
    "\n",
    "  Args:\n",
    "    num_points (int): number of points desired in random sample\n",
    "\n",
    "  Returns:\n",
    "    dataX, dataY (ndarray, ndarray): arrays of size (num_points,) containing x\n",
    "    and y coordinates of sampled points\n",
    "\n",
    "  \"\"\"\n",
    "\n",
    "  # Generate desired number of points uniformly between 0 and 1 (using uniform) for\n",
    "  #     both x and y\n",
    "  dataX = np.random.uniform(0, 1, size = (num_points,))\n",
    "  dataY = np.random.uniform(0, 1, size = (num_points,))\n",
    "\n",
    "  return dataX, dataY\n",
    "\n",
    "# Set a seed\n",
    "np.random.seed(0)\n",
    "\n",
    "# Set number of points to draw\n",
    "num_points = 10\n",
    "\n",
    "# Draw random points\n",
    "dataX, dataY = generate_random_sample(num_points)\n",
    "\n",
    "# Visualize\n",
    "with plt.xkcd():\n",
    "  plot_random_sample(dataX, dataY, \"Random sample of 10 points\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Demo 1.1: Random Sample Generation from Uniform Distribution\n",
    "In practice this may not look very uniform, although that is of course part of the randomness! Uniform randomness does not mean smoothly uniform. When we have very little data it can be hard to see the distribution.\n",
    "\n",
    "Below, you can adjust the number of points sampled with a slider. Does it look more uniform now? Try increasingly large numbers of sampled points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:28.055308Z",
     "iopub.status.busy": "2021-06-02T20:55:28.054728Z",
     "iopub.status.idle": "2021-06-02T20:55:28.323235Z",
     "shell.execute_reply": "2021-06-02T20:55:28.322728Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "\n",
    "#@markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "def generate_random_sample(num_points):\n",
    "  \"\"\" Generate a random sample containing a desired number of points (num_points)\n",
    "  in the range [0, 1] using a random number generator object.\n",
    "\n",
    "  Args:\n",
    "    num_points (int): number of points desired in random sample\n",
    "\n",
    "  Returns:\n",
    "    dataX, dataY (ndarray, ndarray): arrays of size (num_points,) containing x\n",
    "    and y coordinates of sampled points\n",
    "\n",
    "  \"\"\"\n",
    "\n",
    "  # Generate desired number of points uniformly between 0 and 1 (using uniform) for\n",
    "  #     both x and y\n",
    "  dataX = np.random.uniform(0, 1, size = (num_points,))\n",
    "  dataY = np.random.uniform(0, 1, size = (num_points,))\n",
    "\n",
    "  return dataX, dataY\n",
    "\n",
    "@widgets.interact\n",
    "def gen_and_plot_random_sample(num_points = widgets.SelectionSlider(options=[(\"%g\"%i,i) for i in np.arange(0, 2000, 10)])):\n",
    "\n",
    "  dataX, dataY = generate_random_sample(num_points)\n",
    "  fig, ax = plt.subplots()\n",
    "  ax.set_xlabel('x')\n",
    "  ax.set_ylabel('y')\n",
    "  plt.xlim([-0.25, 1.25])\n",
    "  plt.ylim([-0.25, 1.25])\n",
    "  plt.scatter(dataX, dataY)\n",
    "  fig.suptitle(\"Random sample of \" + str(num_points) + \" points\", size=16)\n",
    "  plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 1.2: Random walk\n",
    "\n",
    "Stochastic models can be used to create models of behaviour. As an example, imagine that a rat is placed inside a novel environment, a box. We could try and model its exploration behaviour by assuming that for each time step it takes a random uniformly sampled step in any direction (simultaneous random step in x direction and random step in y direction)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coding Exercise 1.2: Modeling a random walk\n",
    "\n",
    "Use the `generate_random_sample` function from above to obtain the random steps the rat takes at each time step and complete the generate_random_walk function below. For plotting, the box will be represented graphically as the unit square enclosed by the points (0, 0) and (1, 1)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:28.333711Z",
     "iopub.status.busy": "2021-06-02T20:55:28.333118Z",
     "iopub.status.idle": "2021-06-02T20:55:28.342434Z",
     "shell.execute_reply": "2021-06-02T20:55:28.341824Z"
    }
   },
   "outputs": [],
   "source": [
    "def generate_random_walk(num_steps, step_size):\n",
    "  \"\"\" Generate the points of a random walk within a 1 X 1 box.\n",
    "\n",
    "  Args:\n",
    "    num_steps (int): number of steps in the random walk\n",
    "    step_size (float): how much each random step size is weighted\n",
    "\n",
    "  Returns:\n",
    "    x, y (ndarray, ndarray): the (x, y) locations reached at each time step of the walk\n",
    "\n",
    "  \"\"\"\n",
    "  x = np.zeros(num_steps + 1)\n",
    "  y = np.zeros(num_steps + 1)\n",
    "\n",
    "  ###################################################################\n",
    "  ## TODO for students: Collect random step values with function from before\n",
    "  ## Fill out the following then remove\n",
    "  raise NotImplementedError(\"Student exercise: need to complete generate_random_walk\")\n",
    "  ###################################################################\n",
    "\n",
    "  # Generate the uniformly random x, y steps for the walk\n",
    "  random_x_steps, random_y_steps = ...\n",
    "\n",
    "  # Take steps according to the randomly sampled steps above\n",
    "  for step in range(num_steps):\n",
    "\n",
    "    # take a random step in x and y. We remove 0.5 to make it centered around 0\n",
    "    x[step + 1] = x[step] + (random_x_steps[step] - 0.5)*step_size\n",
    "    y[step + 1] = y[step] + (random_y_steps[step] - 0.5)*step_size\n",
    "\n",
    "    # restrict to be within the 1 x 1 unit box\n",
    "    x[step + 1]= min(max(x[step + 1], 0), 1)\n",
    "    y[step + 1]= min(max(y[step + 1], 0), 1)\n",
    "\n",
    "  return x, y\n",
    "\n",
    "# Set a random seed\n",
    "np.random.seed(2)\n",
    "\n",
    "# Select parameters\n",
    "num_steps = 100   # number of steps in random walk\n",
    "step_size = 0.5   # size of each step\n",
    "\n",
    "# Generate the random walk\n",
    "x, y = generate_random_walk(num_steps, step_size)\n",
    "\n",
    "# Visualize\n",
    "plot_random_walk(x, y, \"Rat's location throughout random walk\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:28.366240Z",
     "iopub.status.busy": "2021-06-02T20:55:28.364070Z",
     "iopub.status.idle": "2021-06-02T20:55:28.643176Z",
     "shell.execute_reply": "2021-06-02T20:55:28.642689Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def generate_random_walk(num_steps, step_size):\n",
    "  \"\"\" Generate the points of a random walk within a 1 X 1 box.\n",
    "\n",
    "  Args:\n",
    "    num_steps (int): number of steps in the random walk\n",
    "    step_size (float): how much each random step size is weighted\n",
    "\n",
    "  Returns:\n",
    "    x, y (ndarray, ndarray): the (x, y) locations reached at each time step of the walk\n",
    "\n",
    "  \"\"\"\n",
    "  x = np.zeros(num_steps + 1)\n",
    "  y = np.zeros(num_steps + 1)\n",
    "\n",
    "  # Generate the uniformly random x, y steps for the walk\n",
    "  random_x_steps, random_y_steps = generate_random_sample(num_steps)\n",
    "\n",
    "  # Take steps according to the randomly sampled steps above\n",
    "  for step in range(num_steps):\n",
    "\n",
    "    # take a random step in x and y. We remove 0.5 to make it centered around 0\n",
    "    x[step + 1] = x[step] + (random_x_steps[step] - 0.5)*step_size\n",
    "    y[step + 1] = y[step] + (random_y_steps[step] - 0.5)*step_size\n",
    "\n",
    "    # restrict to be within the 1 x 1 unit box\n",
    "    x[step + 1]= min(max(x[step + 1], 0), 1)\n",
    "    y[step + 1]= min(max(y[step + 1], 0), 1)\n",
    "\n",
    "  return x, y\n",
    "\n",
    "# Set a random seed\n",
    "np.random.seed(2)\n",
    "\n",
    "# Select parameters\n",
    "num_steps = 100   # number of steps in random walk\n",
    "step_size = 0.5   # size of each step\n",
    "\n",
    "# Generate the random walk\n",
    "x, y = generate_random_walk(num_steps, step_size)\n",
    "\n",
    "# Visualize\n",
    "with plt.xkcd():\n",
    "  plot_random_walk(x, y, \"Rat's location throughout random walk\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We put a little green dot for the starting point and a red point for the ending point. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Demo 1.2: Varying parameters of a random walk\n",
    "In the interactive demo below, you can examine random walks with different numbers of steps or step sizes, using the sliders. Make sure to hit \"Run Simulation\" after each change.\n",
    "\n",
    "\n",
    "1.  What could an increased step size mean for the actual rat's movement we are simulating?\n",
    "2. For a given number of steps, is the rat be more likely to visit all general areas of the arena with a big step size or small step size?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:28.691774Z",
     "iopub.status.busy": "2021-06-02T20:55:28.691190Z",
     "iopub.status.idle": "2021-06-02T20:55:28.707394Z",
     "shell.execute_reply": "2021-06-02T20:55:28.700574Z"
    }
   },
   "outputs": [],
   "source": [
    "#@markdown Make sure you execute this cell to enable the widget!\n",
    "def gen_and_plot_random_walk(num_steps, step_size):\n",
    "  x, y = generate_random_walk(num_steps, step_size)\n",
    "  plot_random_walk(x, y, \"Rat's location throughout random walk\")\n",
    "\n",
    "\n",
    "widget=interactive(gen_and_plot_random_walk, {'manual': True},\n",
    "                num_steps=(0, 1000, 1),\n",
    "                step_size=(0, 1, 0.1))\n",
    "\n",
    "widget.children[-2].description='Run Simulation'\n",
    "widget.children[-2].style.button_color='lightgreen'\n",
    "controls = HBox(widget.children[:-1], layout=Layout(flex_flow='row wrap'))\n",
    "output = widget.children[-1]\n",
    "display(VBox([controls, output]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:28.711372Z",
     "iopub.status.busy": "2021-06-02T20:55:28.710790Z",
     "iopub.status.idle": "2021-06-02T20:55:28.713215Z",
     "shell.execute_reply": "2021-06-02T20:55:28.713687Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "1) A larger step size could mean that the rat is moving faster, or that we sample\n",
    "   the rats location less often.\n",
    "\n",
    "2) The rat tends to visit more of the arena with a large step size.\n",
    "\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice a uniform random movement is too simple an assumption. Rats do not move completely randomly; even if you could assume that, you would need to approximate with a more complex probability distribution.\n",
    "\n",
    "Nevertheless, this example highlights how you can use sampling to approximate behaviour. \n",
    "\n",
    "**Main course preview:** On day W3D2 we will see how random walk models can be used to also model accumulation of information in decision making."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 2: Discrete distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:28.718882Z",
     "iopub.status.busy": "2021-06-02T20:55:28.718307Z",
     "iopub.status.idle": "2021-06-02T20:55:28.790927Z",
     "shell.execute_reply": "2021-06-02T20:55:28.791370Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 2: Discrete distributions\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"-_5DtYgwtaA\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The uniform distribution is very simple, and can only be used in some rare cases. If we only had access to this distribution, our statistical toolbox would be very empty. Thankfully we do have some more advanced distributions!\n",
    "\n",
    "The uniform distribution that we looked at above is an example of a continuous distribution. The value of $X$ that we draw from this distribution can take **any value** between $a$ and $b$.\n",
    "\n",
    "However, sometimes we want to be able to look at discrete events. Imagine that the rat from before is now placed in a T-maze, with food placed at the end of both arms. Initially, we would expect the rat to be choosing randomly between the two arms, but after learning it should choose more consistently.\n",
    "\n",
    "A simple way to model such random behaviour is with a single **Bernoulli trial**, that has two outcomes, {$Left, Right$}, with probability $P(Left)=p$ and $P(Right)=1-p$ as the two mutually exclusive possibilities (whether the rat goes down the left or right arm of the maze).\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2.1: Binomial distributions\n",
    "\n",
    "The binomial distribution simulates $n$ number of binary events, such as the $Left, Right$ choices of the random rat in the T-maze. Imagine that you have done an experiment and found that your rat turned left in 7 out of 10 trials. What is the probability of the rat turning left 7 times ($k = 7$)?\n",
    "\n",
    "This is given by:\n",
    "\n",
    "$$ P(k|n,p)= \\left( \\begin{array} \\\\n \\\\ k\\end{array} \\right) p^k (1-p)^{n-k}$$\n",
    "$$\\binom {n}{k}={\\frac {n!}{k!(n-k)!}}$$\n",
    "\n",
    "In this formula, $p$ is the probability of turning left, $n$ is the number of binary events, or trials, and $k$ is the number of times the rat turned left. The term $\\binom {n}{k}$ is the binomial coefficient.\n",
    "\n",
    "If we assume an equal chance of turning left or right, then $p=0.5$. Note that if we only have a single trial $n=1$ this is equivalent to a single Bernoulli trial (feel free to do the math!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Think! 2.1: Binomial distribution sampling\n",
    "We will draw a desired number of random samples from a binomial distribution, with $n = 10$ and $p = 0.5$. Each sample returns the number of trials, $k$, a rat turns left out of $n$ trials.\n",
    "\n",
    "We will draw 1000 samples of this (so it is as if we are observing 10 trials of the rat, 1000 different times). We can do this using numpy: `np.random.binomial(n, p, size = (n_samples,))`\n",
    "\n",
    "See below to visualize a histogram of the different values of $k$, or the number of times the rat turned left in each of the 1000 samples and discuss the following questions.\n",
    "\n",
    "\n",
    "1.   What are the x-axis limits of the histogram and why?\n",
    "2.   What is the shape of the histogram?\n",
    "3.  Looking at the histogram, how would you interpret the outcome of the simulation if you didn't know what p was? Would you have guessed p = 0.5?\n",
    "3.   What do you think the histogram would look like if the probability of turning left is 0.8 ($p = 0.8$)?\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:28.815463Z",
     "iopub.status.busy": "2021-06-02T20:55:28.810534Z",
     "iopub.status.idle": "2021-06-02T20:55:28.983878Z",
     "shell.execute_reply": "2021-06-02T20:55:28.982988Z"
    }
   },
   "outputs": [],
   "source": [
    "# @markdown Execute this cell to see visualization\n",
    "\n",
    "# Select parameters for conducting binomial trials\n",
    "n = 10\n",
    "p = 0.5\n",
    "n_samples = 1000\n",
    "\n",
    "# Set random seed\n",
    "np.random.seed(1)\n",
    "\n",
    "# Now draw 1000 samples by calling the function again\n",
    "left_turn_samples_1000 = np.random.binomial(n, p, size = (n_samples,))\n",
    "\n",
    "# Visualize\n",
    "count, bins = plot_hist(left_turn_samples_1000, 'Number of left turns in sample')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the limits for this histogram? How would you describe the shape?\n",
    "\n",
    "Looking at the histogram, how would you interepret the outcome of the simulation? Do you think it is likely that the rat chose with equal probability?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:28.988408Z",
     "iopub.status.busy": "2021-06-02T20:55:28.987802Z",
     "iopub.status.idle": "2021-06-02T20:55:28.991791Z",
     "shell.execute_reply": "2021-06-02T20:55:28.991026Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "1) The limits of the histogram at 0 and 10, as these are the minimum and maximum\n",
    "number of trials that a rat can turn left out of 10 trials.\n",
    "\n",
    "2) The shape seems symmetric and is centered around 5.\n",
    "\n",
    "3) An average/mean around 5 left turns out of 10 trials indicates that the rat\n",
    "chose left and right turns in the maze with equal probability (that p = 0.5)\n",
    "\n",
    "4) With p = 0.8, the center of the histogram would be at x = 8 and it would not be\n",
    "   as symmetrical (since it would be cut off at max 10). You can go into the code above\n",
    "   and run it with p = 0.8 to see this.\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Note that in the more general case where there are $n$ possible outcomes (our rat is an n-armed maze) each with their own associated probability $p_1, p_2, p_3, p_4, ...$ , we do have to make sure that the probabilities sum to one:\n",
    "\n",
    "$$\\sum_i P(x=i)=\\sum_i p_i =1$$\n",
    "\n",
    "This is the **categorical distribution**; draws from this distribution are a simple extension of the Bernoulli trial.\n",
    "\n",
    "If we sample from this distribution multiple times, we can then describe the distribution of outcomes from each sample as the **multinomial distribution** (generalisation from the binomial distribution)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 2.2: Poisson Distribution\n",
    "For some phenomena there may not be a natural limit on the maximum number of possible events or outcomes. \n",
    "\n",
    "The Poisson distribution is a '**point-process**', meaning that it determines the number of discrete 'point', or binary, events that happen within a fixed space or time, allowing for the occurence of a potentially infinite number of events. The Poisson distribution is specified by a single parameter $\\lambda$ that encapsulates the mean number of events that can occur in a single time or space interval (there will be more on this concept of the 'mean' later!). The formula for a Poisson distribution on x is: \n",
    "\n",
    "$$P(x)=\\frac{\\lambda^x e^{-\\lambda}}{x!}$$\n",
    "\n",
    "Relevant to us, we can model the number of times a neuron spikes within a time interval using a Poisson distribution. In fact, neuroscientists often do! As an example, if we are recording from a neuron that tends to fire at an average rate of 4 spikes per second, then the Poisson distribution specifies the distribution of recorded spikes over one second, where $\\lambda=4$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coding Exercise 2.2: Poisson distribution sampling\n",
    "\n",
    "In the exercise below we will draw some samples from the Poisson distribution and see what the histogram looks.\n",
    "\n",
    "In the code, fill in the missing line so we draw 5 samples from a Poisson distribution with $\\lambda = 4$. Use `np.random.poisson`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:28.996981Z",
     "iopub.status.busy": "2021-06-02T20:55:28.996407Z",
     "iopub.status.idle": "2021-06-02T20:55:28.999687Z",
     "shell.execute_reply": "2021-06-02T20:55:28.999223Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set random seed\n",
    "np.random.seed(0)\n",
    "\n",
    "# Draw 5 samples from a Poisson distribution with lambda = 4\n",
    "sampled_spike_counts = ...\n",
    "\n",
    "# Print the counts\n",
    "print(\"The samples drawn from the Poisson distribution are \" +\n",
    "          str(sampled_spike_counts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.005381Z",
     "iopub.status.busy": "2021-06-02T20:55:29.004805Z",
     "iopub.status.idle": "2021-06-02T20:55:29.008609Z",
     "shell.execute_reply": "2021-06-02T20:55:29.009087Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "\n",
    "# Set random seed\n",
    "np.random.seed(0)\n",
    "\n",
    "# Draw 5 samples from a Poisson distribution with lambda = 4\n",
    "sampled_spike_counts = np.random.poisson(4, 5)\n",
    "\n",
    "# Print the counts\n",
    "print(\"The samples drawn from the Poisson distribution are \" +\n",
    "          str(sampled_spike_counts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should see that the neuron spiked 6 times, 7 times, 1 time, 8 times, and 4 times in 5 different intervals."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Demo 2.2: Varying parameters of Poisson distribution\n",
    "\n",
    "Use the interactive demo below to vary $\\lambda$ and the number of samples, and then visualize the resulting histogram.\n",
    "\n",
    "\n",
    "1.  What effect does increasing the number of samples have?  \n",
    "2.  What effect does changing $\\lambda$ have?\n",
    "3.  With a small lambda, why is the distribution asymmetric?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.039348Z",
     "iopub.status.busy": "2021-06-02T20:55:29.038757Z",
     "iopub.status.idle": "2021-06-02T20:55:29.244581Z",
     "shell.execute_reply": "2021-06-02T20:55:29.243956Z"
    }
   },
   "outputs": [],
   "source": [
    "#@markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "@widgets.interact(lambda_value = widgets.FloatSlider(value=4, min=0, max=10, step=0.1),\n",
    "                  n_samples = widgets.IntSlider(value=5, min=0, max=1000, step=1))\n",
    "\n",
    "def gen_and_plot_possion_samples(lambda_value, n_samples):\n",
    "  sampled_spike_counts = np.random.poisson(lambda_value, n_samples)\n",
    "  count, bins = plot_hist(sampled_spike_counts, 'Recorded spikes per second')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.248225Z",
     "iopub.status.busy": "2021-06-02T20:55:29.247643Z",
     "iopub.status.idle": "2021-06-02T20:55:29.252394Z",
     "shell.execute_reply": "2021-06-02T20:55:29.251897Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove explanation\n",
    "\n",
    "\"\"\"\n",
    "1) More samples means the distribution starts to look more symmetric.\n",
    " This is a result of the Central Limit Theorem, where a large number\n",
    "of sampled values makes samples from the Poisson distribution look similarly\n",
    "distributed to the symmetric Gaussian (Normal) distribution.\n",
    "\n",
    "2) Increasing lambda moves the distribution along the x-axis, essentialy changing the\n",
    "   mean of the distribution. With lambda = 6 for example, we would expect to often see\n",
    "   6 spike counts per interval.\n",
    "\n",
    "3) It becomes asymmetric because the Poisson distribution is restricted to non-negative counts\n",
    "   (you can't have negative numbers of spikes)\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Section 3: Continuous distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.256665Z",
     "iopub.status.busy": "2021-06-02T20:55:29.255786Z",
     "iopub.status.idle": "2021-06-02T20:55:29.334218Z",
     "shell.execute_reply": "2021-06-02T20:55:29.334725Z"
    }
   },
   "outputs": [],
   "source": [
    "#@title Video 3: Continuous distributions\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"nAYbAXDF9h4\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtu.be/\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We do not have to restrict ourselves to only probabilistic models of discrete events. While some events in neuroscience are discrete (e.g. number of spikes by a neuron), many others are continuous (e.g. neuroimaging signals in EEG or fMRI, distance traveled by an animal, human pointing in a direction of a stimulus).\n",
    "\n",
    "With continuous distributions we have to replace the normalising sum $$\\sum_i P(x=p_i) =1$$ over all possible events, with an integral\n",
    "$$\\int_a^b P(x) =1$$ where a and b are the limits of the random variable $X$ (often $-\\infty$ and $\\infty$).\n",
    "\n",
    "If we want to make predictions about possible outcomes (\"I believe the BOLD signal from the area will be in the range $x_1$ to $ x_2 $\") we can use the integral $\\int_{x_1}^{x_2} P(x) =1$.\n",
    "$P(x)$ is now a probability density function.\n",
    "\n",
    "While for discrete outcomes we can ask about the probability of an specific event (\"what is the probability this neuron will fire 4 times in the next second\"), this is not defined for a continuous distribution (\"what is the probability of the BOLD signal being exactly 4.000000000...\"). Hence we need to focus on intervals when calculating probabilities from a continuous distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 3.1: Gaussian Distribution\n",
    "\n",
    "(***Exercise moved  from NMA2020 Bayes day, all credit to original creators!)***\n",
    "\n",
    "The most widely used continuous distribution is probably the Gaussian (also known as Normal) distribution. It is extremely common across all kinds of statistical analyses. Because of the central limit theorem, many quantities are Gaussian distributed. Gaussians also have some nice mathematical properties that permit simple closed-form solutions to several important problems. \n",
    "\n",
    "As a working example, imagine that a human participant is asked to point in the direction where they perceived a sound coming from. As an approximation, we can assume that the variability in the direction/orientation they point towards is Gaussian distributed.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coding Exercise 3.1A: Gaussian Distribution\n",
    "\n",
    "In this exercise, you will implement a Gaussian by filling in the missing portions of code for the function `my_gaussian` below. Gaussians have two parameters. The **mean** $\\mu$, which sets the location of its center, and its \"scale\" or spread is controlled by its **standard deviation** $\\sigma$, or **variance** $\\sigma^2$ (i.e. the square of standard deviation). **Be careful not to use one when the other is required.**\n",
    "\n",
    "The equation for a Gaussian is:\n",
    "$$\n",
    "\\mathcal{N}(\\mu,\\sigma^2) = \\frac{1}{\\sqrt{2\\pi\\sigma^2}}\\exp\\left(\\frac{-(x-\\mu)^2}{2\\sigma^2}\\right)\n",
    "$$\n",
    "Also, don't forget that this is a probability distribution and should therefore sum to one. While this happens \"automatically\" when integrated from $-\\infty$ to $\\infty$, your Gaussian will only be computed over a finite number of points (for the cell below we will sample from -8 to 9 in step sizes of 0.1). You therefore need to explicitly normalize it to sum to one yourself. \n",
    "\n",
    "Test out your implementation with a $\\mu = -1$ and $\\sigma = 1$. After you have it working, play with the parameters to develop an intuition for how changing $\\mu$ and $\\sigma$ alter the shape of the Gaussian. You may want to revisit this when learning about Bayesian statistics in the main course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.349162Z",
     "iopub.status.busy": "2021-06-02T20:55:29.348544Z",
     "iopub.status.idle": "2021-06-02T20:55:29.351286Z",
     "shell.execute_reply": "2021-06-02T20:55:29.351792Z"
    }
   },
   "outputs": [],
   "source": [
    "def my_gaussian(x_points, mu, sigma):\n",
    "  \"\"\" Returns normalized Gaussian estimated at points `x_points`, with\n",
    "  parameters: mean `mu` and standard deviation `sigma`\n",
    "\n",
    "  Args:\n",
    "      x_points (ndarray of floats): points at which the gaussian is evaluated\n",
    "      mu (scalar): mean of the Gaussian\n",
    "      sigma (scalar): standard deviation of the gaussian\n",
    "\n",
    "  Returns:\n",
    "      (numpy array of floats) : normalized Gaussian evaluated at `x`\n",
    "  \"\"\"\n",
    "\n",
    "  ###################################################################\n",
    "  ## TODO for students: Implement the formula for a Gaussian\n",
    "  ## Add code to calculate the gaussian px as a function of mu and sigma,\n",
    "  ## for every x in x_points\n",
    "  ## Function Hints: exp -> np.exp()\n",
    "  ##                 power -> z**2\n",
    "  ##\n",
    "  ## Fill out the following then remove\n",
    "  raise NotImplementedError(\"Student exercise: need to implement Gaussian\")\n",
    "  ###################################################################\n",
    "  px = ...\n",
    "\n",
    "  # as we are doing numerical integration we may have to remember to normalise\n",
    "  # taking into account the stepsize (0.1)\n",
    "  px = px/(0.1*sum(px))\n",
    "  return px\n",
    "\n",
    "x = np.arange(-8, 9, 0.1)\n",
    "\n",
    "# Generate Gaussian\n",
    "px = my_gaussian(x, -1, 1)\n",
    "\n",
    "# Visualize\n",
    "my_plot_single(x, px)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.374356Z",
     "iopub.status.busy": "2021-06-02T20:55:29.358074Z",
     "iopub.status.idle": "2021-06-02T20:55:29.638772Z",
     "shell.execute_reply": "2021-06-02T20:55:29.638267Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "def my_gaussian(x_points, mu, sigma):\n",
    "  \"\"\" Returns normalized Gaussian estimated at points `x_points`, with\n",
    "  parameters: mean `mu` and standard deviation `sigma`\n",
    "\n",
    "  Args:\n",
    "      x_points (ndarray of floats): points at which the gaussian is evaluated\n",
    "      mu (scalar): mean of the Gaussian\n",
    "      sigma (scalar): standard deviation of the gaussian\n",
    "\n",
    "  Returns:\n",
    "      (numpy array of floats) : normalized Gaussian evaluated at `x`\n",
    "  \"\"\"\n",
    "\n",
    "  px = 1/(2*np.pi*sigma**2)**1/2 *np.exp(-(x_points-mu)**2/(2*sigma**2))\n",
    "\n",
    "  # as we are doing numerical integration we may have to remember to normalise\n",
    "  # taking into account the stepsize (0.1)\n",
    "  px = px/(0.1*sum(px))\n",
    "  return px\n",
    "\n",
    "x = np.arange(-8, 9, 0.1)\n",
    "\n",
    "# Generate Gaussian\n",
    "px = my_gaussian(x, -1, 1)\n",
    "\n",
    "# Visualize\n",
    "with plt.xkcd():\n",
    "  my_plot_single(x, px)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Interactive Demo 3.1: Sampling from a Gaussian distribution\n",
    "\n",
    "Now that we have gained a bit of intuition about the shape of the Gaussian, let's imagine that a human participant is asked to point in the direction of a sound source, which we then measure in horizontal degrees. To simulate that we draw samples from a Normal distribution:\n",
    "\n",
    "$$x \\sim \\mathcal{N}(\\mu,\\sigma) $$\n",
    "\n",
    "\n",
    "We can sample from a Gaussian with mean mu and standard deviation sigma using `np.random.normal(mu, sigma, size = (n_samples,))`.\n",
    "\n",
    "In the demo below, you can change the mean and standard deviation of the Gaussian, and the number of samples, we can compare the histogram of the samples to the true analytical distribution (in red).\n",
    "\n",
    "\n",
    "\n",
    "1.   With what number of samples would you say that the full distribution (in red) is well approximated by the histogram? \n",
    "2.   What if you just wanted to approximate the variables that defined the distribution, i.e. mean and variance?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.677804Z",
     "iopub.status.busy": "2021-06-02T20:55:29.676896Z",
     "iopub.status.idle": "2021-06-02T20:55:29.919054Z",
     "shell.execute_reply": "2021-06-02T20:55:29.918460Z"
    }
   },
   "outputs": [],
   "source": [
    "#@markdown Make sure you execute this cell to enable the widget!\n",
    "\n",
    "\n",
    "@widgets.interact(mean = widgets.FloatSlider(value=0, min=-5, max=5, step=0.5),\n",
    "                  standard_dev = widgets.FloatSlider(value=0.5, min=0, max=10, step=0.1),\n",
    "                  n_samples = widgets.IntSlider(value=5, min=0, max=1000, step=1))\n",
    "def gen_and_plot_normal_samples(mean, standard_dev, n_samples):\n",
    "  x = np.random.normal(mean, standard_dev, size = (n_samples,))\n",
    "  xspace = np.linspace(-10, 10, 100)\n",
    "  plot_gaussian_samples_true(x, xspace, mean, standard_dev,\n",
    "                            'orientation (degrees)', 'probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Main course preview:** Gaussian distriutions are everywhere and are critical for filtering, linear systems (W2D2), optimal control (W3D3) and almost any statistical model of continuous data (W3D1, W3D2, etc.)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "Across the different exercises you should now:\n",
    "* have gotten some intuition about how stochastic randomly generated data can be\n",
    "* understand how to model data using simple distributions \n",
    "* understand the difference between discrete and continuous distributions\n",
    "* be able to plot a Gaussian distribution\n",
    "\n",
    "For more reading on these topics see:\n",
    "Textbook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# Bonus: Markov chains\n",
    "\n",
    "If you have the time, here is an extra exercise.\n",
    "\n",
    "Next we examine a example Markov chain.\n",
    "Imagine that a rat is able to move freely between 3 areas: a dark rest area\n",
    "(1), a nesting area (2) and a bright area for collecting food (3). Every 5 minutes we record the rat's location. The table below shows the probability of the rat transitioning from one area to another ($state_i$ to $state_{i+1}$).\n",
    "\n",
    "\\begin{array}{|l |  l | l | l |} \\hline\n",
    "state_{i} &P(state_{i+1}=1|state_i=*) &P(state_{i+1}=2|state_i=*) & P(state_{i+1}=3|state=_i*) \\\\ \\hline\n",
    "state_{i}=1& 0.2  &0.6 &0.2\\\\\n",
    "state_{i}=2& .6 &0.3& 0.1\\\\\n",
    "state_{i}=3& 0.8 &0.2 &0\\\\ \\hline\n",
    "\\end{array}\n",
    "\n",
    "We are modeling this as a Markov chain, so the animal is only in one of the states at a time and can transition between the states.\n",
    "\n",
    "If the animal starts in area 2, what is the probability the animal will again be in area 2 when we check on it 20 minutes (4 transitions) later?\n",
    "\n",
    "Fill in the transition matrix below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.929638Z",
     "iopub.status.busy": "2021-06-02T20:55:29.927528Z",
     "iopub.status.idle": "2021-06-02T20:55:29.932472Z",
     "shell.execute_reply": "2021-06-02T20:55:29.931856Z"
    }
   },
   "outputs": [],
   "source": [
    "###################################################################\n",
    "## TODO for student: Fill transition matrix with table values\n",
    "## Fill out the following then remove\n",
    "raise NotImplementedError(\"Student exercise: need to fill transition matrix\")\n",
    "###################################################################\n",
    "\n",
    "# Fill out transition matrix with matching table values\n",
    "# in the form np.array([ [...], [...], [...] ])\n",
    "transition_matrix = ...\n",
    "\n",
    "# Initial state, x0\n",
    "x0 = np.array([0, 1, 0])\n",
    "\n",
    "# Then multiply initial state by the matrix 4 times, or equivalently raise the\n",
    "# matrix to the 4th power before multiplying it to the initial state\n",
    "x4 = x0 @ np.linalg.matrix_power(transition_matrix, 4)\n",
    "\n",
    "# The second area is indexed as 1 (Python starts indexing at 0)\n",
    "print(\"The probability the rat will be in area 2 after 4 transitions is: \" + str(x4[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.938880Z",
     "iopub.status.busy": "2021-06-02T20:55:29.937617Z",
     "iopub.status.idle": "2021-06-02T20:55:29.940790Z",
     "shell.execute_reply": "2021-06-02T20:55:29.940293Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "\n",
    "# Fill out transition matrix with matching table values\n",
    "# in the form np.array([ [...], [...], [...] ])\n",
    "transition_matrix = np.array([[ 0.2, 0.6, 0.2],[ .6, 0.3, 0.1], [0.8, 0.2, 0]])\n",
    "\n",
    "# Initial state, x0\n",
    "x0 = np.array([0, 1, 0])\n",
    "\n",
    "# Then multiply initial state by the matrix 4 times, or equivalently raise the\n",
    "# matrix to the 4th power before multiplying it to the initial state\n",
    "x4 = x0 @ np.linalg.matrix_power(transition_matrix, 4)\n",
    "\n",
    "# The second area is indexed as 1 (Python starts indexing at 0)\n",
    "print(\"The probability the rat will be in area 2 after 4 transitions is: \" + str(x4[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should get a probability of 0.4311, i.e. there is a 43.11% chance that you will find the rat in area 2 in 20 minutes.\n",
    "\n",
    "What is the average amount of time spent by the rat in each of the states Implicit in the question is the idea that we can start off with a random initial state and then measure how much relative time is spent in each area. If we make a few assumptions (e.g. ergodic system) we can instead start with an initial random distribution and see how that changes after 100 time steps/transitions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.949546Z",
     "iopub.status.busy": "2021-06-02T20:55:29.948962Z",
     "iopub.status.idle": "2021-06-02T20:55:29.952190Z",
     "shell.execute_reply": "2021-06-02T20:55:29.951665Z"
    }
   },
   "outputs": [],
   "source": [
    "# Initialize random initial distribution\n",
    "x_random = np.ones((1,3))/3\n",
    "\n",
    "###################################################################\n",
    "## TODO for student: Fill compute the state matrix after 100 transitions\n",
    "raise NotImplementedError(\"Student exercise: need to complete computation below\")\n",
    "###################################################################\n",
    "\n",
    "# Fill in the missing line to get the state matrix after 100 transitions, like above\n",
    "x_average_time_spent = ...\n",
    "print(\"The proportion of time spend by the rat in each of the three states is: \"\n",
    "            + str(x_average_time_spent[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-02T20:55:29.958241Z",
     "iopub.status.busy": "2021-06-02T20:55:29.956923Z",
     "iopub.status.idle": "2021-06-02T20:55:29.960413Z",
     "shell.execute_reply": "2021-06-02T20:55:29.959813Z"
    }
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "\n",
    "# Initialize random initial distribution\n",
    "x_random = np.ones((1,3))/3\n",
    "\n",
    "# Fill in the missing line to get the state matrix after 100 transitions, like above\n",
    "x_average_time_spent = x_random @ np.linalg.matrix_power(transition_matrix, 100)\n",
    "print(\"The proportion of time spend by the rat in each of the three states is: \"\n",
    "            + str(x_average_time_spent[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The proportion of time spend in each of the three areas are 0.4473, 0.4211, and 0.1316, respectively.\n",
    "\n",
    "If the animal is satiated and tired the transitions change to:\n",
    "\n",
    "\\begin{array}{|l |  l | l | l |} \\hline\n",
    "state_{i} &P(state_{i+1}=1|state_i=*) &P(state_{i+1}=2|state_i=*) &P(state_{i+1}=3|state_i=*) \\\\ \\hline\n",
    "state_{i}=1& 0.2  &0.7 &0.1\\\\\n",
    "state_{i}=2& .3 &0.7& 0.\\\\\n",
    "state_{i}=3& 0.8 &0.2 &0\\\\ \\hline\n",
    "\\end{array}\n",
    "\n",
    "Try repeating the questions above for this table of transitions by changing the transition matrix. See how much time the rat spends on average in each area.\n",
    "What would you predict?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Main course preview:** The Markov property is extremely important for many models, particularly Hidden Markov Models, discussed on day W3D2, and for methods such as Markov Chain Monte Carlo sampling."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W0D5_Tutorial1",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
